Intro to BioEnergeticFoodWebs

by Chris Griffiths, Eva Delmas and Andrew Beckerman, Dec. 2020.

7.7 μs

This document follows on from 'Getting started', 'Basic Julia commands' and 'Differential Equations in Julia' and assumes that you're still working in your active project.

This document introduces the BioEnergeticFoodWebs.jl and EcologicalNetworks.jl packages. It demonstrates how to run the BioEnergetic Food Web (BEFW) model, how to vary variables of interest (e.g., productivity) and construct experiments designed to investigate the effect of different variables on population and community dynamics.

For those that are unfamilar with the BEFW and it's application in Julia, we advise checking out the MEE paper before we start. Remember, the BEFW model is also based on a system of differential equations and is solved using the same engine as the DifferentialEquations.jl package.

7.0 μs

Load packages

3.0 μs

You'll need the following packages for this tutorial:

2.7 μs

xxxxxxxxxx
 
using BioEnergeticFoodWebs, EcologicalNetworks, JLD2, Statistics, Plots, CSV, DataFrames, Random

36.4 s

The JLD2.jl package will be useful later as it allows you to directly export and load a BEFW output object. Let's also set a random seed for reproducibility:

3.1 μs

MersenneTwisterseedUInt321

0x00000015

stateRandom.DSFMT.DSFMT_statevalInt321

163718196

1073255743

-1451494363

1072899844

-1618160924

1073258805

-1517321307

1073339839

-1150166203

1072942963

-1583246675

1072950051

-407139083

1073171397

1426608413

1073390650

507415766

1073502426

-1219197164

1073024304

more761

-1705970718

762

1072820091

763

394465429

764

1073380663

765

887692544

766

1617232906

767

-1051299616

768

1513956225

769

770

valsFloat641

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

0.0

more993

0.0

994

0.0

995

0.0

996

0.0

997

0.0

998

0.0

999

0.0

1000

0.0

1001

0.0

1002

0.0

intsUInt1281

0x00000000000000000000000000000000

0x00000000000000000000000000000000

0x00000000000000000000000000000000

0x00000000000000000000000000000000

0x00000000000000000000000000000000

0x00000000000000000000000000000000

0x00000000000000000000000000000000

0x00000000000000000000000000000000

0x00000000000000000000000000000000

0x00000000000000000000000000000000

0x00000000000000000000000000000000

0x00000000000000000000000000000000

0x00000000000000000000000000000000

0x00000000000000000000000000000000

0x00000000000000000000000000000000

0x00000000000000000000000000000000

0x00000000000000000000000000000000

0x00000000000000000000000000000000

0x00000000000000000000000000000000

0x00000000000000000000000000000000

more492

0x00000000000000000000000000000000

493

0x00000000000000000000000000000000

494

0x00000000000000000000000000000000

495

0x00000000000000000000000000000000

496

0x00000000000000000000000000000000

497

0x00000000000000000000000000000000

498

0x00000000000000000000000000000000

499

0x00000000000000000000000000000000

500

0x00000000000000000000000000000000

501

0x00000000000000000000000000000000

idxF

idxI

xxxxxxxxxx
 
Random.seed!(21)

122 ms

Preamble

One of main advantages of running food web models in Julia is that simulations are fast and can be readily stored in your active project. With this in mind, make a new folder in your project called out_objects (right click > New Folder). Alternatively, you can create an out_objects folder directly using mkdir().

4.9 μs

"example_folder"

xxxxxxxxxx
 
# We've already create a folder called out_objects in our project but an example of mkdir() would be:
mkdir("example_folder")
# if you haven't created an out_objects folder yet, simply replace "example_folder" with "out_objects". 

101 μs

Running the BEFW

There are four major steps when running the BioEnergetic Food Web model in Julia:

Generate an initial network
Fix parameters
Simulate
Explore output and plot

Initial network

Before running the BEFW model, we have to construct an initial random network using the niche model. The network is characterised by the number of species in the network and its connectance value. Here, we generate a network of 20 species with a connectance value of 0.15:

11.6 μs

20×20 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0
 1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮              ⋮           
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0

x
 
begin
    # generate network
    A_bool = EcologicalNetworks.nichemodel(20,0.15) 
    # convert the UnipartiteNetwork object into a matrix of 1s and 0s
    A = Int.(A_bool.A)
end

36.2 μs

In the above code chunk, we are saving the output from running the nichemodel as A_bool and then using the A part of A_bool to construct our initial random network. Within A, 1s indicate an interaction among species and 0s no interaction. In the packages used here, the networks are directed from i (rows) to j (columns), describing the direction of the interaction (i eats j), not of the flow of biomass.

8.1 μs

You can check the connectance of A using:

2.8 μs

0.1975

xxxxxxxxxx
 
# calculate connectance
co = sum(A)/(size(A,1)^2)

1.2 μs

Here, connectance is calculated as the number of realised links (sum of 1s in A) divided by the number of species in A squared. This end part (species^2) describes the maximum number of possible links in the network A.

7.1 μs

Parameters

Prior to running the BEFW model, you have to create a vector of model parameters using the model_parameters function. Numerous parameter values can be specified within the model_parameters function, however, most of them have default values that are built into the BioEnergeticFoodWebs.jl package. For simplicity, we use the default values here:

7.4 μs

Dict

:α

1.0

:e_carnivore

0.85

:Γh

Float641

0.0

0.0

0.0

0.5

0.0

0.5

0.5

0.5

0.5

0.5

0.5

0.5

0.5

0.5

0.0

0.5

0.5

0.5

0.5

0.5

:m_producer

1.0

:c

0.0

:h

1.0

:vertebrates

BitArray{1}: 1

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

:tmpA

Any

:rewire_method

:none

:trophic_rank

Float641

1.0

1.0

1.0

2.0

1.0

2.71429

2.25

2.75

2.25

2.66667

2.33333

3.0

3.5

2.85714

1.0

3.1

2.66667

3.1

2.66667

3.0

:w

20×20 Array{Float64,2}:
 0.0   0.0   0.0       0.0       0.0        …  0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   1.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.142857  0.142857  0.142857   …  0.0       0.0       0.0  0.0  0.0
 0.25  0.25  0.25      0.25      0.0           0.0       0.0       0.0  0.0  0.0
 ⋮                                          ⋱  ⋮                             
 0.0   0.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.1        …  0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.166667  0.166667  0.166667      0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.1           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0           0.333333  0.333333  0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0909091     0.0       0.0       0.0  0.0  0.0

:TSR_type

:no_response

:e_herbivore

0.45

:productivity

:species

:efficiency

20×20 Array{Float64,2}:
 0.0   0.0   0.0   0.0   0.0   0.0   0.0   …  0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.45  0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.45  0.85  0.45  0.85  0.85  …  0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.45  0.45  0.45  0.85  0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 ⋮                             ⋮           ⋱              ⋮                     
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.45  0.85  0.85  …  0.85  0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.45  0.85  0.45  0.85  0.85     0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.45  0.85  0.85     0.85  0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.45  0.85  0.85  0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.45  0.85  0.85     0.85  0.45  0.0   0.0   0.0  0.0  0.0

:K

1.0

:S

:extinctionstime

Any

:is_producer

BitArray{1}: 1

true

true

true

false

true

false

false

false

false

false

false

false

false

false

true

false

false

false

false

false

:dp

#10

:x

Float641

0.138

0.138

0.138

0.3141

0.138

0.3141

0.3141

0.3141

0.3141

0.3141

0.3141

0.3141

0.3141

0.3141

0.138

0.3141

0.3141

0.3141

0.3141

0.3141

:Z

1.0

:extinctions

Any

:dry_mass_293

Float641

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

:np

:A

20×20 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0
 1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮              ⋮           
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0

:Γ

Float641

0.0

0.0

0.0

0.5

0.0

0.5

0.5

0.5

0.5

0.5

0.5

0.5

0.5

0.5

0.0

0.5

0.5

0.5

0.5

0.5

:ar

Float641

0.0

0.0

0.0

16.0

0.0

16.0

16.0

16.0

16.0

16.0

16.0

16.0

16.0

16.0

0.0

16.0

16.0

16.0

16.0

16.0

:bodymass

Float641

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

:is_herbivore

BitArray{1}: 1

false

false

false

true

false

true

true

true

true

true

true

false

false

true

false

true

true

true

true

true

xxxxxxxxxx
 
# create model parameters
p = model_parameters(A)
# in the most simple case, the model_parameters function simply requires A

746 μs

For more information and a full list of the parameters and their defaults values type ?model_parameters in the REPL.

If you want to view, check or extract any of the parameter values in p use the p[:name] notation. For example, you can view a vector of each species' trophic rank using:

7.1 μs

Float641

1.0

1.0

1.0

2.0

1.0

2.71429

2.25

2.75

2.25

2.66667

2.33333

3.0

3.5

2.85714

1.0

3.1

2.66667

3.1

2.66667

3.0

xxxxxxxxxx
 
# view trophic ranks:
p[:trophic_rank]

439 ns

Simulate

To run the BEFW model, we first assign biomasses at random to each species and then simulate the biomass dynamics forward using the simulate function:

5.0 μs

Dict

:p

Dict

:α

1.0

:e_carnivore

0.85

:Γh

Float641

0.0

0.0

0.0

0.5

0.0

0.5

0.5

0.5

0.5

0.5

0.5

0.5

0.5

0.5

0.0

0.5

0.5

0.5

0.5

0.5

:m_producer

1.0

:c

0.0

:h

1.0

:vertebrates

BitArray{1}: 1

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

:tmpA

Any1

20×20 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0
 1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮              ⋮           
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0

20×20 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0
 1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮              ⋮           
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0

20×20 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0
 1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮              ⋮           
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0

20×20 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0
 1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮              ⋮           
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0

20×20 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0
 1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮              ⋮           
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0

20×20 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0
 1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮              ⋮           
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0

20×20 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0
 1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮              ⋮           
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0

20×20 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0
 1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮              ⋮           
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0

20×20 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0
 1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮              ⋮           
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0

20×20 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0
 1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮              ⋮           
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0

:rewire_method

:none

:trophic_rank

Float641

1.0

1.0

1.0

2.0

1.0

2.71429

2.25

2.75

2.25

2.66667

2.33333

3.0

3.5

2.85714

1.0

3.1

2.66667

3.1

2.66667

3.0

:w

20×20 Array{Float64,2}:
 0.0   0.0   0.0       0.0       0.0        …  0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   1.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.142857  0.142857  0.142857   …  0.0       0.0       0.0  0.0  0.0
 0.25  0.25  0.25      0.25      0.0           0.0       0.0       0.0  0.0  0.0
 ⋮                                          ⋱  ⋮                             
 0.0   0.0   0.0       0.0       0.0           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.1        …  0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.166667  0.166667  0.166667      0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.1           0.0       0.0       0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0           0.333333  0.333333  0.0  0.0  0.0
 0.0   0.0   0.0       0.0       0.0909091     0.0       0.0       0.0  0.0  0.0

:TSR_type

:no_response

:e_herbivore

0.45

:productivity

:species

:efficiency

20×20 Array{Float64,2}:
 0.0   0.0   0.0   0.0   0.0   0.0   0.0   …  0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.45  0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.45  0.85  0.45  0.85  0.85  …  0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.45  0.45  0.45  0.85  0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 ⋮                             ⋮           ⋱              ⋮                     
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.45  0.85  0.85  …  0.85  0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.45  0.85  0.45  0.85  0.85     0.0   0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.45  0.85  0.85     0.85  0.0   0.0   0.0   0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.45  0.85  0.85  0.0  0.0  0.0
 0.0   0.0   0.0   0.0   0.45  0.85  0.85     0.85  0.45  0.0   0.0   0.0  0.0  0.0

:K

1.0

:S

:extinctionstime

Any11

22.2043

29.861

55.4288

75.6713

103.35

178.837

322.529

598.673

94.5363

101

152.679

:is_producer

BitArray{1}: 1

true

true

true

false

true

false

false

false

false

false

false

false

false

false

true

false

false

false

false

false

:dp

#10

:x

Float641

0.138

0.138

0.138

0.3141

0.138

0.3141

0.3141

0.3141

0.3141

0.3141

0.3141

0.3141

0.3141

0.3141

0.138

0.3141

0.3141

0.3141

0.3141

0.3141

:Z

1.0

:extinctions

Any1

:dry_mass_293

Float641

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

0.153846

:np

:A

20×20 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0
 1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮              ⋮           
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  1  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  0  0  0
 0  0  0  0  1  1  1  1  1  1  1  1  1  1  1  0  0  0  0  0

:Γ

Float641

0.0

0.0

0.0

0.5

0.0

0.5

0.5

0.5

0.5

0.5

0.5

0.5

0.5

0.5

0.0

0.5

0.5

0.5

0.5

0.5

:ar

Float641

0.0

0.0

0.0

16.0

0.0

16.0

16.0

16.0

16.0

16.0

16.0

16.0

16.0

16.0

0.0

16.0

16.0

16.0

16.0

16.0

:bodymass

Float641

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

:is_herbivore

BitArray{1}: 1

false

false

false

true

false

true

true

true

true

true

true

false

false

true

false

true

true

true

true

true

:B

8021×20 Adjoint{Float64,Array{Float64,2}}:
 0.034165   0.175658    0.677504  0.571024   …  0.50232   0.928288   0.0639268
 0.036847   0.0451099   0.434938  0.433516      0.634629  1.12844    0.0804242
 0.0423065  0.0149496   0.318599  0.322718      0.776915  1.32988    0.0978262
 0.0506449  0.00659167  0.26418   0.248923      0.9215    1.4966     0.115019
 0.0617965  0.00361401  0.241252  0.199829      1.05942   1.59266    0.130852
 0.0759696  0.00232627  0.236711  0.165936   …  1.18146   1.60765    0.144321
 0.0935981  0.00168823  0.24454   0.141479      1.27958   1.56103    0.154667
 ⋮                                           ⋱                       
 0.299553   0.200912    0.258188  0.0126103  …  0.0       0.0311834  0.79528
 0.298484   0.200092    0.25727   0.0126249     0.0       0.0311948  0.795905
 0.29731    0.199179    0.256248  0.0126282     0.0       0.0312078  0.796613
 0.296065   0.1982      0.255152  0.0126202     0.0       0.0312223  0.797385
 0.294786   0.197183    0.254016  0.0126011     0.0       0.0312377  0.798198
 0.293516   0.196166    0.252877  0.0125716  …  0.0       0.0312535  0.799025

:t

Float641

0.0

0.25

0.5

0.75

1.0

1.25

1.5

1.75

2.0

2.25

2.5

2.75

3.0

3.25

3.5

3.75

4.0

4.25

4.5

4.75

more8012

1997.75

8013

1998.0

8014

1998.25

8015

1998.5

8016

1998.75

8017

1999.0

8018

1999.25

8019

1999.5

8020

1999.75

8021

2000.0

xxxxxxxxxx
 
begin
    # assign biomasses
    bm = rand(size(A,1)) 
    # select biomasses at random between ]0:1[
    # as an alternative, you could assign all species the same biomass of 1 using bm = ones(size(A,1)) 
    
    # simulate
    out = simulate(p, bm, start=0, stop=2000)
    # this might take a few seconds
end

390 ms

The simulate function requires the model parameters p and the species biomasses bm. In addition, you can specify the timespan of the simulation (using the start and stop arguments), fix a species extinction threshold (using extinction_threshold) and select a solver (using use). For more information type ?simulate in the REPL.

4.9 μs

Output and plot

Once the simulation finishes, the output is stored as a dictionary called out. Within out there are three entries:

out[:p] - lists the parameters
out[:B] - biomass of each species through time
out[:t] - timesteps (these typically increase in 0.25 intervals)

The biomass dynamics of each species can then be plotted. Similar to the DifferentialEquations.jl package, the BioEnergeticFoodWebs.jl package also has it's own built in plotting recipe:

11.3 μs

xxxxxxxxxx
 
# plot
Plots.plot(out[:t], out[:B], legend = true, ylabel = "Biomass", xlabel = "Time")
# this may take a minute to render

10.1 ms

You'll notice that the biomass dynamics are noisey during the first few hundred time steps, these are the system's transient dynamics. The dynamics then settle into a steady state where the system can be assumed to be at equilbirum. You'll also notice that some species go extinct and some persist, the initial number of species in the food web (20 in this case) can found using out[:p][:S] and the identity of those that went extinct using out[:p][:extinctions].

The BioEnergeticFoodWebs.jl package also has a range of built in functions that conveniently calculate some of the key metrics of the food web, these include the total biomass, the diversity, the species persistence and the temporal stability:

7.7 μs

biomass

2.398808427844779

xxxxxxxxxx
 
# total biomass
biomass = total_biomass(out, last=1000)

62.2 μs

diversity

0.8310935582838079

xxxxxxxxxx
 
# diversity
diversity = foodweb_evenness(out, last=1000)

450 μs

persistence

0.5

xxxxxxxxxx
 
# persistence
persistence = species_persistence(out, last=1000)

249 μs

stability

-0.012570492028439493

xxxxxxxxxx
 
# stability 
stability = population_stability(out, last=1000)

166 μs

Each of these functions will output a single value. This value is the average over the last 1000 time steps. For more information, use ? to access the help files on each function in the REPL (e.g., ?species_persistence).

3.5 μs

Variables

Once you've got the BEFW model running, the next step is to vary a variable of interest and rerun. For example, we might be interested in what affect a small change in Z (consumer-resource body mass ratio) has on the estimated food web and its biomass dynamics. The default value for Z is 1.0, but what happens if we increase it to 10.0:

4.5 μs

x
 
begin
    # set Z (has to be a floating number not an integer)
    Z = 10.0
    # create model parameters
    p_z = model_parameters(A, Z = Z)
    # assign biomasses
    bm_z = rand(size(A,1)) 
    # simulate
    out_z = simulate(p_z, bm_z, start=0, stop=2000)
    # plot
    Plots.plot(out_z[:t], out_z[:B], legend = true, ylabel = "Biomass", xlabel = "Time")
end

213 ms

Similarly, what happens if we also increase the carrying capacity (K) of the resource from 1.0 (default) to 5.0:

3.1 μs

xxxxxxxxxx
 
begin
    # set K (has to be a floating number not an integer)
    K = 5.0 
    # create model parameters
    p_K = model_parameters(A, Z = Z, K = K)
    # assign biomasses
    bm_K = rand(size(A,1)) 
    # simulate
    out_K = simulate(p_K, bm_K, start=0, stop=2000)
    # plot
    Plots.plot(out_K[:t], out_K[:B], legend = true, ylabel = "Biomass", xlabel = "Time")
end

555 ms

As you've probably guessed, the main message here is that many variables can be changed in the BEFW model and it's super easy to do so. Some changes will have large effects and some not so much. In the next step, we take this one step further.

6.5 μs

Experiments

The next step is to construct a computional experiment designed to investigate the effect of different variables on population and community dynamics. To do this we construct a gradient of variables as vectors and then simulate the BEFW model multiple times using a loop. To illustrate this, we're going to reproduce example 1 from Delmas et al. 2016. The aim of this example is to investigate the effect of increasing K on food web diversity. In addition, we're also going to allow α (interspecific competition relative to intraspecific competition) to vary and repeat the experiment 5 times with 5 different initial networks.

First, we define the experiment by creating vectors of our variables and fixing the number of repetitions:

5.8 μs

Float641

0.92

1.0

1.08

xxxxxxxxxx
 
α = [0.92, 1.0, 1.08]
# 0.92 - the interspecific competition is smaller than the intraspecific competition promoting coexistence
# 1.0 - neutrally stable 
# 1.08 - the intraspecific competition is smaller the interspecific competition favouring competitive exclusion

645 ns

Float641

0.1

0.16681

0.278256

0.464159

0.774264

1.29155

2.15443

3.59381

5.99484

10.0

xxxxxxxxxx
 
# vector of K
k = exp10.(range(-1,1,length=10))
# log scale from 0.1 to 10

49.6 ms

reps

xxxxxxxxxx
 
# number of reps
reps = 5

35.0 ns

We then create a dataframe to store the outputs:

2.9 μs

	α	K	network	diversity	stability	biomass
	Any	Any	Any	Any	Any	Any

xxxxxxxxxx
 
# dataframe
df = DataFrame(α = [], K = [], network = [], diversity = [], stability = [], biomass = [])

167 ms

and construct a while loop to generate the 5 unique initial networks, each of which contains 20 species with a connectance value of 0.15:

3.2 μs

xxxxxxxxxx
 
begin
    # list to store networks
    global networks = []
    # monitoring variable 
    global l = length(networks)
    # while loop
    while l < reps
        # generate network
        A_bool = EcologicalNetworks.nichemodel(20,0.15) 
        # convert the UnipartiteNetwork object into a matrix of 1s and 0s
        A = Int.(A_bool.A)
        # calculate connectance
        co = sum(A)/(size(A,1)^2)
        # ensure that connectance = 0.15
        if co == 0.15
            push!(networks, A)
            # save network is co = 0.15
        end
        global l = length(networks)
    end
end

1.9 ms

We can then run the simulations by looping, using nested for loops, over the unique values of α and K, as well as the 5 unique initial networks. After each simulation we will save each output object to our active project as a JLD2 file and store any output metrics of interest in our dataframe:

3.3 μs

xxxxxxxxxx
 
# loop over networks
for h in 1:reps
    A = networks[h]
    # here, you might want to save a copy of the initial network using writedlm(A)
​
    # loop over α
    for i in 1:length(α) 
        # loop over K
        for j in 1:length(k)
​
        # create model parameters
​
        p = model_parameters(A, α = α[i], K = k[j])
        # assign biomasses
        bm = rand(size(A,1)) 
        # simulate
        out = simulate(p, bm, start=0, stop=2000)
​
        # dummy naming variables
        α_num = α[i]
        K_num = k[j]
        # save `out` as a JLD2 object using the @save macro:
        @save "out_objects/model_output, network = $h, alpha = $α_num, K = $K_num.jld2" out
​
        # calculate output metrics
        diversity = foodweb_evenness(out, last = 1000)
        stability = population_stability(out, last = 1000)
        biomass = total_biomass(out, last = 1000)
​
        # push to df
        push!(df, [α[i], k[j], h, diversity, stability, biomass])
​
        # print some stuff - see how the simulation is progressing
        println(("α = $α_num", "K = $K_num", "network = $h"))
        end
    end
end
# the code will be much faster if you remove the @save command

98.6 s

We can then explore the outputs and plot our results. Here, instead of using the built in plotting recipe, we will construct a plot that matches figure 1 in Delmas et al. 2016. Specifically, we will plot food web diversity (y-axis) as a function of K (x-axis) and α (colour):

4.0 μs

	variable	mean	min	median	max	nunique	nmissing	eltype
	Symbol	Float64	Float64	Float64	Float64	Int64	Int64	DataType
1	:α	1.0	0.92	1.0	1.08	3	0	Any
2	:K	2.48181	0.1	1.03291	10.0	10	0	Any
3	:network	3.0	1.0	3.0	5.0	5	0	Any
4	:diversity	0.833252	0.336135	0.843676	1.0	150	0	Any
5	:stability	-0.180435	-1.50381	-0.0165022	-7.14755e-16	150	0	Any
6	:biomass	2.88808	0.301331	2.12444	12.8342	150	0	Any

xxxxxxxxxx
 
# explore output
describe(df)

1.0 s

	α	K	network	diversity	stability	biomass
	Any	Any	Any	Any	Any	Any
1	0.92	0.1	1.0	0.932747	-0.000169469	0.301331
2	0.92	0.16681	1.0	0.938526	-0.000168663	0.426691
3	0.92	0.278256	1.0	0.83412	-0.00014034	0.600914
4	0.92	0.464159	1.0	0.825871	-0.000139997	0.970888
5	0.92	0.774264	1.0	0.82981	-0.27333	1.63807
6	0.92	1.29155	1.0	0.964897	-0.000777308	1.05171

xxxxxxxxxx
 
first(df,6)

17.0 μs

	α	K	network	diversity	stability	biomass
	Any	Any	Any	Any	Any	Any
1	1.08	0.774264	5.0	0.846921	-0.102008	2.71432
2	1.08	1.29155	5.0	0.880867	-0.358162	3.58376
3	1.08	2.15443	5.0	0.842566	-0.606752	4.05996
4	1.08	3.59381	5.0	0.889399	-0.407231	4.05905
5	1.08	5.99484	5.0	0.649198	-1.16421	5.54389
6	1.08	10.0	5.0	0.999372	-0.00794146	1.8543

xxxxxxxxxx
 
last(df,6)

14.8 μs

xxxxxxxxxx
 
# plot
# initialise an empty plot
pl = Plots.plot([NaN], [NaN],
                label = "",
                ylims = (0,1.1),
                leg = :bottomright,
                foreground_colour_legend = nothing,
                xticks = (log10.(k), string.(round.(k, digits = 1))),
                xlabel = "Carrying capacity",
                ylabel = "Food web diversity (evenness)")

345 ms

shp

Symbol1

:square

:diamond

:utriangle

xxxxxxxxxx
 
# set marker shapes
shp = [:square, :diamond, :utriangle]

900 ns

Symbol1

:solid

:dash

:dot

xxxxxxxxxx
 
# set line types
ls = [:solid, :dash, :dot]

652 ns

clr

xxxxxxxxxx
 
# set colours
clr = [RGB(174/255, 139/255, 194/255), RGB(188/255, 188/255, 188/255), RGB(124/255, 189/255, 122/255)]
# when we define colours in Julia they are printed 

18.2 ms

lbl

String1

"Coexistence"

"Neutral"

"Exclusion"

xxxxxxxxxx
 
# set legend labels
lbl = ["Coexistence", "Neutral", "Exclusion"]

867 ns

xxxxxxxxxx
 
# make the plot
for (i, α) in enumerate(α)
    # subset
    tmp = df[df.α .== α, :]
    # remove NaN values
    tmp = tmp[.!(isnan.(tmp.diversity)), :]
    # calculate mean across reps
    meandf = by(tmp, :K, :diversity => mean)
    # command to avoid printing legends multiple times
    l = i == 1 ? lbl[i] : ""
    # add to pl
    plot!(pl, log10.(meandf.K), meandf.diversity_mean,
              msc = clr[i],
              mc = :white,
              msw = 3,
              markershape = shp[i],
              linestyle = ls[i],
              lc = clr[i],
              lw = 2,
              label = lbl[i],
              seriestype = [:line :scatter])
end

3.6 s

xxxxxxxxxx
 
# display plot
plot(pl)

56.2 ms

Finally, we can save our dataframe as a .csv file:

2.9 μs

"My_data.csv"

xxxxxxxxxx
 
# save
CSV.write("My_data.csv", df)

682 ms